Submitted to Phys. Rev. A 

Photon Statistics; Nonlinear Spectroscopy of Single Quantum 

Systems 

Shaul Mukamel 

Department of Chemistry, University of California, Irvine, CA 92697-2025 

(Dated: February 2, 2008) 

Abstract 

A unified description of multitime correlation functions, nonlinear response functions, and quan- 
tum measurements is developed using a common generating function which allows a direct com- 
parison of their information content. A general formal expression for photon counting statistics 
from single quantum objects is derived in terms of Liouville space correlation functions of the ma- 
terial system by making a single assumption that spontaneous emission is described by a master 
equation. 
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I. INTRODUCTION 



Knowledge of the full density matrix of the radiation field allows to compute all its 
measurable properties. In particular, photon counting statistics which had proven to be 
a most valuable measure of coherence has been formulated in the sixties by Kelley and 
Kleiner, Glauber, and Mandel in terms of expectation values of normally-ordered field op- 
erator a 1 i 2 i?i 4 i^^7i§ . A revived interest had emerged in the eighties when stochastic trajectory 
experiments on single quantum objects, atoms, ion traps, molecules and quantum dots be- 
came feasibl o 9ll0lllll2ll3i:L4115116117 . It is more convenient to formulate such measurements in 
terms of correlation functions of the material system, rather than of the radiation field. 
Various applications for specific few-level model systems have been investigated 9 - ' 18 ' 19 ' 20 ' 21 . 

In this paper we develop correlation function expressions for photon statistics which ap- 
ply to a general model of a quantum system driven by an external field and coupled to a 
bath. The normally-ordered field expressions are remarkably general; the only assumption 
made in their derivation is that the photon detection is described by the Fermi golden rule. 
Similarly, our material correlation function expressions hold under a single assumption; that 
spontaneous emission can be described by a master equation. The derivation of the master 
equation starting with the fully quantum description of the field is well documente d 22 ' 23 ' 24 . 
No other properties of the radiation field enter explicitly in the present formulation. Com- 
puting the reduced density matrix of a single quantum system coupled to a bath has been a 
long standing goal of nonequilibrium statistical mechanic a 24 ' 25 ' 26 . Various types of reduced 
equations of motion based on stochastic or microscopic models are well developed. The 
present approach is therefore particularly useful for single quantum systems since it can uti- 
lize any level of reduced equations of motion to predict the photon statistics. Computing the 
many-body density matrix of a macroscopic system is much more complex and a collective 
description using field operators^' 3 ^'^'^ may then be more adequate. 

We further develop some fundamental connections between multitime correlation func- 
tions of photon statistics and response functions of nonlinear spectroscopy^. Coherent 
experiments conducted using multiple pulses provide a wealth of information on electronic 
and nuclear dynamics 28 . These techniques can create and manipulate quantum coherences 
among selected states and the signals provide snapshots of their dynamics. There are nu- 
merous motivations for performing such measurements, (i) Novel spectroscopy: the ability 
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to explore and access unusual regions of phase space, (ii) Coherent control: achieving a de- 
sired goal (e.g. optimize the branching ratio of a reaction towards a favorable product), (hi) 
Quantum computing: these applications depend on generating and retrieving information 
about coherences between several degrees of freedom prepared in correlated wavepackets 
(such entanglement is a synonym to the old fashioned term correlation), (iv) Overcome 
coupling to a bath e.g. selectively eliminating dephasing processes(or the more trendy term 
decoherence) . We shall draw upon the analogy between photon statistics and nonlinear 
response and correlation functions to show important similarities and differences in their 
information content and simulation strategies. 

In Section II we present the Liouville space expressions for multipoint correlation func- 
tions and response functions. In Section III we discuss multipoint quantum equilibrium 
measurements and introduce a unified generating function that can be used to compute 
correlation, response and measurements. This sets the stage for deriving the Liouville space 
expressions for photon counting in Section IV. Finally, our results are summarized in Section 
V. 

II. MULTIPOINT GENERATING FUNCTIONS FOR CORRELATION AND RE- 
SPONSE FUNCTIONS 

The state of complex quantum systems may be conveniently characterized by multitime 
quantities which carry various levels of information and are easier to calculate, measure, 
or visualize compared to the many-body wavefunction^. In this section we present formal 
expressions for two such objects (correlation and response functions) using a Liouville space 
(superoperator) approach. While the results of this section are not new, they establish the 
notation, setting the stage for calculating the successive measurements in Section III and 
eventually to the photon statistics in Section IV, which are our main results. Correlation 
and response functions can be defined in Hilbert space using ordinary operators. Quantum 
measurements on the other hand, must be formulated using density matrices in Liouville 
space. The notation introduced here is essential for expressing all of these quantities us- 
ing a common generating function (Eq. (j2SI)) which unambiguously reveals their relative 
information content. 

Consider a dynamical variable of interest A which can represent e.g. the dipole operator, 
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the coordinate or the momentum of a tagged particle, or some collective coordinate. The 
simplest n-point object is the correlation function 

C {n \r n . . . r„) = Tr[A(r n ) . . . A{ n ) Peq l (1) 

where p eq is the equilibrium density matrix of the system. 

Classical correlation functions are given by moments of the joint distribution of successive 
measurements and are therefore directly observable. Quantum correlation functions, in 
contrast, are not connected to specific measurements in a simple way. Instead, response 
function o 25 i 27 i 29 i 30 which represent the reaction of the system to an external field E(t) coupled 
to the variable A via Hi nt = —E(t) ■ A may be readily measured. In a response experiment 
the total Hamiltonian Ht{t) consists of material Hamiltonian H and the coupling to the 
driving field 

H T (T) = H + H mt (r)- (2) 
We shall be interested in the expectation value of A at time t 

U(t) = Tr[Ap(t)} = ((A\p(t)}), (3) 

where \p(t))) = J2jk Pikif)\jk)) is the density matrix of the system, and the ket \ jk)) denotes 
the Liouville-space operator \j)(k\ 25 i 27 i 31 . 
Eq. (JHJ) can be recast in the forr^~ ' ' 



U(t) = (TA + (t)exp 



% -f_JrE{r)A^T)\). (4) 



Here and below (• ■ ■) denotes averaging with respect to the equilibrium density matrix p eq 

(A)=Tr[A Peq l (5) 

A± are superoperators acting in Liouville space defined as follows: For any ordinary operator 
A we define 

A _ = A L -A R ; A + = ^ (A L + A R ) , (6) 

where Ai and Ar are the superoperators that act on the ket (left) and bra (right) of the 
density matrix (A L B = AB and A R B = -BA R ). 

The time evolution in Eq. (j3J) is given in the interaction picture. For an ordinary operator 
in Hilbert space this is defined by 

A(t) =exp^#r) Aexpf-^ffr) . (7) 



Similarly, the time evolution of superoperators is governed by the Liouville operator H 
corresponding to the material Hamiltonian H 

Aj(r) = exp (~#-t) A j exp (- jr#-r) . J = +, - L, R (8) 

T is the positive time ordering operator which rearranges all products of superoperators in 
order of decreasing time from left to right. The nonlinear response functions are obtained 
by expanding the exponent of Eq. in powers of E{r). The expectation value of A to 
n — l'th order in the field is given by 

U {n ~ l \r n )= f n dt Tn _, . . . r dT 1 R^(T n ...T 1 )£(r n - 1 )...£(n), n = 2,3... (9) 

J —oo J —oo 

with the nonlinear response function 

R^\r n ... n) ee (£f (A + (r n )A4r n ^) . . . A_(n)). (10) 

ft(n) re p re sents the response at times r n to n — 1 very short pulses applied at times 
Ti ■ ■ -r n _i— . Note that all time arguments are fully ordered T\ < r 2 . . . < r n . The op- 
erator A + {r n ) corresponds to the observation time, the operators A_{rj) j = 1, • • ■ ,n — 1 
represents interactions with the external field at times tj. We chose to label the response 
function corresponding to by R^ rather than R^ n ~^ since it is an n point function; 

this will facilitate the comparison with the other multitime quantities discussed below. 
Eq. (fTUj) is an abbreviated notation for 

R [n \r n ...r 1 )= (£f ([. . . [[A(r n ), A(r n ^)l A(r n _ 2 )] . . . , A( Tl )}), (11) 

or 

R (n) (r n .. .n) = (^j n Tr{A(r n )[A(r n ^), [A fa), [A{r x ),p eq \] ■■•]}. (12) 

R^ is thus given by a combination of n-order ordinary (Hilbert space) correlation func- 
tions. Eq. ()12j) contains 2 n ~ 1 terms representing all possible "left" and "right" actions of 
the various commutators. Each term corresponds to a Liouville- space path and can be rep- 
resented by a double-sided Feynman diagram^. The various pathways interfere, giving rise 
to many interesting effects such as new resonances. For a multilevel system R^ is usually 
expanded in the eigenstates of the free Hamiltonian H. Each path then consists of n — 1 
periods of free evolution separated by n couplings with A, which change the state of the 
system. 



Using our superoperator notation, the correlation function Eq. (JTJ) can be recast as 

C^{r n ...r 1 ) = Tr[A L (r n ) . . . A l {t x )p^ (13) 

where the time evolution of A^ir) is given by Eq. (jSJ). 

Classical quantities are conveniently represented as moments of some joint distribution 
functions connected to measurements. The closest we can come up in quantum mechanics 
is through moments of generating functions. This is not only a convenient computational 
tool, but also provides insights and helps connect different measurements. The generating 
function for correlation functions is defined as 

W^(a n r n ■ ■ ■ am) = (5(a n - A{r n )) ■ ■ ■ 5( ai - A(n))). (14) 

By comparing Eq. (^) and Eq. ()14|) we immediately find that 

C (ri) (r n . . . n) = J ... J ax... a n Wc\a n T n ■ ■ ■ am) da x ... da n . (15) 

Similarly we can define a generating function for response functions 

W^\a n r n ■ ■ ■ am) (16) 

= ([8(a n - A{r n )), ■ ■ ■ [6(a 2 - A(r 2 )), [S( ai - A( n )), Peq }))\ 

so that the response function is given by 

R {n) (r n . . . n) = J . . . J at--- a n W^\a n r n ■ ■ ■ am) da!--- da n . (17) 

Eqs. (fT4^) and (fTo^) play the role of a classical distribution functions even though they are 
generally complex and may be negative. Nevertheless, they serve as generating functions 
for correlation and response functions which are given by their first "moments", Eqs. (fTo^l 
and (fT7f) . 

So far we considered two n point objects: Correlation functions and response functions. 
A third type of quantity which is more closely connected to photon statistics is the joint 
distribution of n successive measurements. This will be introduced next. 

III. UNIFIED GENERATING FUNCTION FOR CORRELATION, RESPONSE 
AND EQUILIBRIUM MEASUREMENTS 

We consider a sequence of n measurements of a dynamical variable A performed on the 
system at times t% < r 2 < r 3 • ■ ■ < r n and yielding the outcomes a x • • • a n . We would like to 
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compute the following ensemble average over many such measurements 

S^r n ■ ■ ■ n = A(r n )---A(r 1 ). (18) 

This quantity is common to classical and quantum systems alike. Eq ()18)1 is a shorthand 
notation for 

S^\r n ---r l )= (19) 

/-/oi-an^W.,-,^) d ai ...da n 

where the joint distribution function Wg 1 ^ can be computed using the theory of quantum 
measurement o? 4 ^^ 1 ^ 1 ^ 1 ^ . To that end, we define the eigenstates \cej > of A with eigenvalues 
a,j, and represent the operator A in the form of an expansion in projection operators Aj 

A = Y,a 3 A J , (20) 

3 

with 

Aj = \aj><a d \. (21) 
We next define the Liouville space projection operator onto the diagonal elements of A 

P(a)=X>(a-a*)l^><^l- (22) 

j 

The Liouville space bracket <C F\G ^>= Tr'FG denotes a scalar product computed by a par- 
tial trace over the measured degrees of freedom of the operator A. Using this projection, the 
joint distribution function of successive measurements may be recast in the f or m 34 i 35 i 36 i 37 i 39 

Ws(a n T n ■ ■ ■ ain) = Tr P(a n , r n ) • • ■ P(oi, n)p eq (23) 

where the time evolution of P(a, r) is given by the interaction picture (Eq. (jHJ)). The compact 
Liouville space notation used in Eq. (}2"3*|) will help establish the connection between photon 
counting and other multitime quantities. Note that both Wc and Ws are normalized as 

J Wj n \a n r n ■ ■ ■ atT\) da\ ■ ■ ■ da n = 1 j — C,S (24) 



whereas Wr, which represents the deviation of the density matrix from equilibrium, has a 

wj^ {a n T n ■ ■ ■ a±Ti) dai ■ ■ ■ da n = 0. (25) 



zero trace 4 ^ 



So far, we introduced three different generating functions to describe correlation func- 
tions, response functions and the joint probability of measurements (Eqs. 1)14)1 . 1)16)1 and 
respectively). To establish the connection between these quantities it will be useful to com- 
pute all three using a single fundamental object. This is possible by using the following 
generating function, 

W( n \a n a' n T n , a 2 a' 2 r 2 ■ ■ ■ a^Ti) = 

Tr {5(a n - A(r n )) ■ ■ ■ 5( ai - A(r 1 ))p e9 5(a' 1 - A( n )) ■ ■ ■ 5(a' n - A(r n ))} . (26) 
The density matrix underlying Eq. ()26|) is 

p [n \a n a' n r n ■ ■ ■ axa'^i) = ^ an< ^ an _ ia ^ 1 (r n - r n _i) • • • Q^a',^ (t 2 - 7i)p aia / (n), (27) 
where 

G(t) = e-xri l — 

h 



^( r )=exp(-li7_T), (28) 



is the interaction-picture propagator. 

The generating function for correlation functions (Eq. (JHJ) is recovered by integrating 
Eq. (|26|) over the primed variables 

Wc\a n T n , ■ ■ ■ , aiTi) = y ^^(ana^rn, a 2 a' 2 r 2 , ■ ■ ■ aia^n) da[da 2 ■ ■ ■ da' n (29) 

and the correlation function is given by Eq. (|15j) . Similarly, the response function is obtained 
by the following integration 

R {n) (r n ■ ■ • n) = J dai-'-ddnJ da[ ■ ■ ■ da' n (a n - a' n ) ■ • ■ (a x - a'^W^ (a n a' n r n ■ ■ ■ a^Ti). 

(30) 

Comparing Eq. ()26|) with Eq. (J2*3*j) . it is clear that the joint probability of measurements 
is related to the diagonal elements of Eq. (}2T)j) . i.e., aj = a'j. However we cannot simply set 
cij = a'j in Eq. (|2*d^) since it will diverge. In order to properly obtain Eq. (126)) from Eq. ()23|) 
we need to add a finite resolution for the measurement defined by a normalized function 
f{a — a') sharply peaked at zero. We can then write 

W^ n \a n r n ---a 1 T 1 ) = (31) 
da[ ■ ■ ■ da'JViata^n, a 2 a' 2 r 2 ■ ■ ■ a n a' n r n )f(a l - a[) ■ ■ ■ f(a n - a' n ) 



Note that the definition of W s is more clean in Liouville space (Eq. (|23j) ) but requires some 
care in Hilbert space. 

We can now better appreciate the fundamental differences between these multitime vari- 
ous quantities. Wc (and C' n ') depend only on cij, where a'j are integrated out. This follows 
from Eq. (fT3|) which only contains "left" superoperators. Because of the — a'- factors 
in Eq. (JHUj) R^ n \ on the other hand, depends only on the off diagonal elements of W with 
a,j 7^ a'j (diagonal elements do not contribute). Finally, Ws (and S^) depends solely on 
the diagonal elements of W with aj = a'j. Ws is the basic quantity in the consistent history 
description of quantum dynamic a?9i?7i?? and has all the properties of a classical joint proba- 
bility distribution. At each time Tj the system is in the state | otj > and its density matrix 
is | a.j >< otj | j = 0, ... n. In contrast, in a nonlinear response measurement as described 
by Ry 1 ' we only measure A at the last time r n ; at the earlier times Tj(j = 1 ... n — 1) we 
only "pass through" otj, but the density matrix could be either | a,j >< ot^ | or | >< ctj \ 
with k 7^ j. Wr is thus not a joint probability; even though we perform some operation on 
the system at n points (n — 1 interactions with the fields and the time of observation), only 
the last interaction corresponds to an actual measurement. In the other times we merely 
perturb the system. It should be emphasized that even though the response functions (and 
W R n ^) are experimental observables that may be obtained from multiple pulse experiments 
with heterodyne detection^., they may not be represented by the joint probability since 
it does not carry enough information for representing this type of observables. 

The response function carries information that depends on delicate interferences among 
events that occur at the various points in time. This interference may be understood in 
terms of sums over pathways which differ by their time ordering i.e., {A{t\)A{t-2)A{t^)), 
(Afa) A(t{) A(t$)) etc. It is less obvious that a similar interference does exist in classical 
mechanics as well. Classically, of course, time ordering is immaterial since all operators 
commute and it suffices to calculate (A{ji) Afa) A{jz)) for n < r 2 • • • < r„. Quantum 
mechanically, each of the n\ permutations of the time arguments in an n point correlation 
function is distinct and carries a different information. Classical correlation functions then 
carry less information than their quantum counterparts. The classical interference takes 
place between closely lying trajectories^*^. 
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IV. CORRELATION FUNCTION EXPRESSIONS FOR PHOTON STATISTICS 



Photon counting statistics as well as shot noise statistics of electrons^ are most closely 
related to (or S^) since they involve n real measurements. However, photon counting 
is a more complex operation than described by since it is performed under nonequi- 
librium conditions where the system is strongly driven by an external field. Furthermore, 
the material system is not closed since photons are emitted. Thirdly, the measurement does 
change the state of the system, not by merely projecting onto a diagonal element. All of 
these complications can be adequately addressed and Eq. (J23j) can be modified to account 
for photon statistics, as will be shown below. 

To proceed further we introduce the master equation, derived by tracing the density 
matrix over the radiation fiel d 22 i 23 i 24 . We shall denote by Y vv i the spontaneous emission rate 
from state v' to the lower energy state v . The total decay rate (inverse radiative lifetime) 
of state v' will then be 

lv , = ]T rv. (32) 

v+v' 

The effects of spontaneous emission are then incorporated by the master equation 

= + 1v)pvv\ v±i/ (33) 

dpw \ -« -p 

—j7~ — —IvPvv + 1 w'Pv'u'- 

Adopting Liouville space (tetradic) notation, the master equation reads 



with 



|=-7P-FP, (34) 



and 

r = J2 \ vv » < v ' v '\- ( 36 ) 

The quantities defined in the previous section correspond to a closed-system and may be 
described either in Hilbert space or in Liouville spac o 44 ' 45 ; the choice is a matter of conve- 
nience. The master equation allows us to avoid the explicit treatment of the field degrees of 
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freedom when describing an open system. This can only be done for the density matrix in 
Liouville space. 

We consider a single quantum system coupled to a bath, driven by an external field and 
subjected to spontaneous emission. The Liouville equation will be partitioned as 



(37) 



Here the Liouville operator L(t) = H-(t) includes our single multilevel system, any other 
bath degrees of freedom, as well as the driving field. It also includes the 7 matrix (Eq. (|35p) 
which represents the diagonal (in Liouville space) part of the master equation. T, on the 
other hand (Eq. (|HEJ)), is the off diagonal part of the master equation which describes the 
transitions among states v and v' . 

We define the Green function solution of Eq. (|H7|) with T = 



</(t 2 ,ti) = Texp 



1 r T2 



drL{r) 



n 



(38) 



and introduce the off diagonal radiative-decay operator in the interaction picture 

r(r) =G\t, 0)r£(r,0). 
The solution of Eq. (j3*Tj) in the interaction picture then reads 



(39) 



p(t) = Texp 



( <1tT{t) 

Jta 



Pit,). 



(40) 



The total probability density of emitting a photon between t and t + dt and another 
photon between t and t + dt is 



W(t,t ) =Tr 



T T(t) exp 



- f drT(r) 

Jt 



r(t )p(*o) 



(41) 



Expanding the solution of Eq. (j4T)j) to n'th order in T yields 



(42) 



where 



n=Q 



P [n ~ l \r n )= r dr n _i... Hdn 
Jo Jo 



G{r n , r„_i)r • • • TQ(t2, ti)Tp(tx). 



(43) 
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pO- 1 ) describes n — 1 photon emission processes at times T\ ■ ■ ■ r n _i. 

The probability density of emitting n photons at times T\ • • • r n is obtained by 

multiplying the integrand by T from the left and taking a trace. This gives 

K {n \r n ■ ■ ■ n) = Tr [T(r ft ) • ■ ■ r(r 1 )p(r 1 )] . (44) 

This general expression for photon statistics in terms of system variables is equivalent to 
the standard normally-ordered expression of field variable a 1 ' 2 ' 3 ^ 6 . To see that, we consider 
the probability of emitting n photons between times ti and tf 

Pn(tf, U) = f' dr n r dr n _i • ■ ■ r dr x K^ n \r n ■■■r l ) (45) 

As can be seen from Eq. (I4UJ) . this is normalized as 

oo 

Y,Pn(t f ,t l )=Trp(t) = 1 (46) 

n=0 

To compute P n we introduce the generating function 

oo 

G(t f ,U;U) = Y,Pn(t f ,U)U n (47) 

n=0 

It then follows from Eqs. (jSJ) and Eq. (gHJ) that 



G(t f ,ti,U) = (Texp 
G thus satisfies the equation of motion. 



U [ tf rfrr(r) 
Ju 



p(U] 



dG{t ^ U) = ~L(t)G(t, U,V)- UTG(t, U, U) 



Using Eq. PSJ) we have 



and Eq. (jlTjl gives 



dW 



-G(U) 



u=o 



= J2Pn(t,t + T)- 

u=x n=0 [n-m)\ 
= (n(n — 1) • • • (n — m + 1)) 



(48) 



(49) 



(50) 



(51) 



which is the n'th factorial moment of P n . Setting m = 1 and m = 2 in Eq. (joTj) we get 

cZG(Z7) 



(n) 



(n 2 ) ~ (n) 



dU 



u=i 



(52) 



d 2 G{U) 
dUo 



u=i 
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Higher moments may be calculated similarly. 

The most commonly used measure of photon statistics, the Mandel parameter, has been 
shown to be related to if t 2 ) 1 ! 9 ? 20 in simple kinetic models of single quantum systems. Photon 
statistics has been calculated using the simplest reduced descriptions such as the Bloch 
equations^. The present approach opens up the use of a broad class of reduced descriptions 
and stochastic models^ for computing K^ 2 \ One interesting application will be to photon 
statistics in superradiance of aggregates^. Eq. (jUj) may be easily generalized to describe 
more refined, frequency resolved, measurements whereby at each time we monitor a different 
(preselected) transition. This can be done simply by using a different element of T at each 
time T(tj) to represent the desired transition. Eq. (JH|) could then provide more detailed 
information about the system. 



V. DISCUSSION 

We have introduced several types of multipoint functions commonly used in experimental 
observations and their theoretical analysis. Using the Liouville space superoperators nota- 
tion, we can recast these various quantities in a formally similar form that facilitates their 
comparison. Eq. (JT3*j) can be written as 

C {n \r n ■ ■ -n) = Tr[A L g(r n - r n _i)A L ■ ■ ■ G(r 2 - r x )A L p^ (53) 

where Q{j) is given by Eq. f)28j) . The nonlinear response function (Eq. (jlOjl ) can be similarly 
recast in the form 

R (n) (r n ■ ■ -n) = Ti[A + g(r n , -r^A-Q ■ • • - n )A- Peq \. (54) 

The joint distribution of successive measurements (Eq. (jUIl)) is written as 

(a n r n ■ ■ ■ ain) = Tr P(a n )g(r n - r n _ 1 )P(a n _ 1 ) • • • P{a 2 )Q{T 2 - T 1 )P{a x )p eq . (55) 

The probability density of observing consecutive photons (Eq. ()44|) ) is 

K {n) (r n • • • n) = Tr [Tg(r n , r„_i) • • • Tg(r 2 , r 1 )r P (r 1 )] (56) 

where g(r, r') is given by Eq. (f3"H|) . Finally, the probability density of measuring n photons 
at times T\ ■ ■ ■ r n (regardless of how many photons are emitted in between) is 

P {n \r n ■ ■ ■ ri) = Tr \rg(r n , r n _0 • • • Qfa, ^Tpfa)] (57) 
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where Q is the Green function solution of Eq. ()H7j) for the driven system 

p(t)=G(t,t )p(t ). (58) 

We note several marked differences between the photon statistics observables (Eqs. (|5fij) 
and (|57jl) and the other quantities (Eq. (|53*)l . and Since the latter are equilibrium 

properties, the Green function is translationally invariant and only depends on the time 
difference Q{jj — r^) rather than on Tj and separately Qir^r^). Also the initial density 
matrix p{j\) in photon statistics measurements is generally a more complex object than p eq 
since it requires computing the preparation stage leading to a nonequilibrium steady state. 
This does not cause any problem in stochastic models where the bath evolution does not 
depend on the state of the system, piji) is then completely specified since the first photon 
emission at T\ determines the state of the system (the final state of the emission) and the 
bath is always in equilibrium. However, fully microscopic modelling will require a separate 
calculation of p[j\). 

Eq. (|5fi|) is very similar to the general expression for n successive measurements (Eq. (|55jl ). 
However, the T matrix is off diagonal since photon emission is accompanied by a transition 
in the system, as opposed to the diagonal P(a) in Eq. (|B3j) which represents ordinary mea- 
surements. Were we to use a diagonal r = \vv ^>^C vv\ it would represent the probability 
of measuring the system at state v at times T\ ■ ■ -r n . Photon counting, however, implies 
that the system is at state v prior to the count but it changes to state v' after the count; 
this is the initial state for the next period of propagation. Apart from this, Eq. (|5fijl or 
Eq. (157)1 are equivalent to n point measurements (Eq. (|55)l). These differences stem from 
the nonequilibrium nature of photon counting performed on open driven systems. 

Finally we note that Eq. ()56)) and Eq. (jSTj) are reminiscent of the normally ordered 
expressions with field operators^ where V represents the detector rather than spontaneous 
emission. In the present approach we do not need normal ordering since in Liouville space 
time ordering is enough to maintain the bookkeeping of interactions. We also note that L(r) 
in Eq. ([38)1 contains the 7 matrix and the Green function therefore contains some diagonal 
signatures of the photon emission. This is required for maintaining the trace of the density 
matrix. Such terms should also be present in the field formulation, but are usually neglected 
and the Green function represents the pure system (without the detector)^. Adding these 
corrections could improve the standard theory of photon statistics. 
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